/** Fichier qui implémente une classe représentant les matrices réelles 3x3.
 *
 * Ce fichier contient l'implémentation de la classe Matrice3x3.
 *
 * @file Matrice3x3.cpp
 * @author Johan "Solo" Luisier
 *
 * @date 24 Juin 2011
 */
#include "Matrice3x3.hpp"

/*
 * Permet de manipuler une ligne / colonne d'une matrice comme un vecteur.
 */
#include <boost/numeric/ublas/matrix_proxy.hpp>

/*
 * Permet d'afficher une matrice simplement dans un flux de sortie.
 */
#include <boost/numeric/ublas/io.hpp>

#include <iostream>

using namespace std;
using namespace boost::numeric::ublas;

/** Définition du symbole de Levi-Civita \f$\varepsilon^{ijk}\f$, le
 * tenseur totalement antisymétrique.
 *
 * J'ai cherché un moyen plus compact de calculer la parité de la permutation
 * créée par les indices, sans succès. Vu le petit nombre de cas à traiter ce
 * n'est pas un problème.
 *
 * @param[in] i premier indice dans le tenseur
 * \f$\varepsilon^{ijk}\f$.
 * @param[in] j deuxième indice dans le tenseur
 * \f$\varepsilon^{ijk}\f$.
 * @param[in] k troisième indice dans le tenseur
 * \f$\varepsilon^{ijk}\f$.
 *
 * @return 1 si \f$\left( i, j, k \right)\f$ est une permutation paire
 * de \f$\left( 1, 2, 3 \right)\f$, -1 pour une permutation impaire et
 * 0 dans tous les autres cas.
 *
 * @author Johan "Solo" Luisier
 *
 * @date 24 Juin 2011
 */
double LeviCivita( const unsigned& i, const unsigned& j, const unsigned k )
{
    if ( ( i == 0 && j == 1 && k == 2 ) ||
         ( i == 1 && j == 2 && k == 0 ) ||
         ( i == 2 && j == 0 && k == 1 ) )
        return 1.;
    else if ( ( k == 0 && j == 1 && i == 2 ) ||
              ( k == 1 && j == 2 && i == 0 ) ||
              ( k == 2 && j == 0 && i == 1 ) )
        return -1.;
    else
        return 0.;
}

/** Constructeur par défaut, construit une matrice scalaire.
 *
 * Une matrice scalaire est égale à la matrice identité multipliée par un
 * nombre quelconque.
 *
 * @param[in] scalaire nombre réel qui se retrouvera sur la diagonale.
 */
Matrice3x3::Matrice3x3( const double& scalaire )
    : matrix< double >( identity_matrix< double > ( 3, 3 ) )
{
    (*this) *= scalaire;
}

/** Construit une matrice réelle à partir de trois vecteurs.
 *
 * Si l'onvoit la matrice comme un endomorphisme de \f$ \mathbb{R} \f$, alors
 * les trois vecteurs sont les images des vecteurs de la base canonique.
 *
 * @param[in] xPrime image de \f$ \vec{e}_1 \f$.
 * @param[in] yPrime image de \f$ \vec{e}_2 \f$.
 * @param[in] zPrime image de \f$ \vec{e}_3 \f$.
 */
Matrice3x3::Matrice3x3( const boost::numeric::ublas::vector< double >& xPrime,
                        const boost::numeric::ublas::vector< double >& yPrime,
                        const boost::numeric::ublas::vector< double >& zPrime )
    : matrix< double >( zero_matrix< double > ( 3, 3 ) )
{
    if ( xPrime.size() != 3 ||
         yPrime.size() != 3 ||
         zPrime.size() != 3 )
    {
        cerr << "Les vecteurs en argument doivent avoir dimension 3 !" << endl;
        exit( 1 );
    }

    matrix_column< matrix< double > > matCol0( *this, 0 );
    matCol0.assign( xPrime );
    matrix_column< matrix< double > > matCol1( *this, 1 );
    matCol1.assign( yPrime );
    matrix_column< matrix< double > > matCol2( *this, 2 );
    matCol2.assign( zPrime );
}

/** Constructeur permettant d'initialiser une matrice diagonale.
 *
 * Une matrice diagonale est telle que @f$ M_{i, j} = 0 \forall i \ne j @f$.
 *
 * @param[in] m11 élément \f$ 1, 1 \f$ de la matrice.
 * @param[in] m22 élément \f$ 2, 2 \f$ de la matrice.
 * @param[in] m33 élément \f$ 3, 3 \f$ de la matrice.
 */
Matrice3x3::Matrice3x3( const double& m11, const double& m22,
                        const double& m33 )
    : matrix< double >( zero_matrix< double > ( 3, 3 ) )
{
    operator()( 0, 0 ) = m11;
    operator()( 1, 1 ) = m22;
    operator()( 2, 2 ) = m33;
}

/** Constructeur de copie.
 *
 * La copie est effectuée par la classe mère.
 *
 * @param[in] mat instance de Matrice3x3 à dupliquer.
 */
Matrice3x3::Matrice3x3( const Matrice3x3& mat )
    : matrix< double >( mat )
{}

/** Destructeur, virtuel pour assurer une destruction propre.
 *
 * Le destructeur n'a rien à faire, la gestion est déléguée à la
 * classe mère.
 */
Matrice3x3::~Matrice3x3()
{}

/** Opérateur d'affectation.
 *
 * Cet opérateur permet de copier une instance dans l'instance courante. Le
 * mécanisme de copie est laissé à la charge de la classe mère.
 *
 * @param[in] mat instance de Matrice3x3 à copier.
 *
 * @return une référence sur l'instance courante.
 */
Matrice3x3& Matrice3x3::operator=( const Matrice3x3& mat )
{
    return operator=( static_cast< MatriceBase >( mat ) );
}

/** Opérateur d'accès à un élément.
 *
 * Il s'agit en fait d'un alias, la méthode matrice< double >::operator() étant
 * déjà surchargée en une version const et une version permettant la
 * modification.
 *
 * @param[in] i indice "ligne".
 * @param[in] j indice "colonne".
 *
 * @return une référence constante sur l'élément \f$ \left( i, j \right) \f$ de
 * la matrice.
 */
const double& Matrice3x3::at( const unsigned& i, const unsigned& j ) const
{
    return operator()( i, j );
}

/** Calcule le déterminant de la matrice.
 *
 * On utilise la définition générale du déterminant : \f[ \det \left(
 * M \right) = \sum_{\sigma} \sum_{i,j,k} \epsilon \left( \sigma
 * \right) \, M_{1, \sigma\left( i \right)} \, M_{2, \sigma\left( j
 * \right)} \, M_{3, \sigma\left( k \right)} \f]
 *
 * @return la valeur du déterminant de la matrice.
 */
double Matrice3x3::determinant() const
{
    double det( 0. );

    unsigned i, j, k;

    for ( i = 0; i < 3; i++ )
    {
        for ( j = 0; j < 3; j++ )
        {
            for ( k = 0; k < 3; k++ )
                det += LeviCivita( i, j, k ) *
                        at( 0, i ) * at( 1, j ) * at( 2, k );
        }
    }

    return det;
}

/** Calcule la trace de la matrice.
 *
 * La trace est simplement la somme des éléments diagonaux.
 *
 * @return la valeur de la trace de la matrice.
 */
double Matrice3x3::trace() const
{
    double resultat( 0. );

    for ( unsigned i( 0u ); i < 3u; i++ )
        resultat += at( i, i );

    return resultat;
}

/** Calcule le cofacteur \f$\left( i, j \right)\f$ de la matrice.
 *
 * Le cofacteur est le déterminant de la sous-matrice \f$\left( i, j
 * \right)\f$.
 *
 * On construit donc la sous-matrice \f$\left( i, j \right)\f$ en
 * annulant la ligne i et la colonne j, puis on assigne \f$ M_{i, j} =
 * 1 \f$. On appelle ensuite la méthode déterminant pour cette matrice
 * intermédiaire.
 *
 * @param[in] i indice "ligne".
 * @param[in] j indice "colonne".
 *
 * @return le cofacteur \f$\left( i, j \right)\f$.
 */
double Matrice3x3::cofacteur( const unsigned& i, const unsigned& j ) const
{
    Matrice3x3 temp( *this );

    VecteurBase col( 3 ), ligne( 3 );

    col = zero_vector< double >( 3 );
    ligne = zero_vector< double >( 3 );

    matrix_column< matrix< double > > matCol( temp, j );
    matrix_row< matrix< double > > matLigne( temp, i );

    matCol.assign( col );
    matLigne.assign( ligne );

    temp( i, j ) = 1.;

    return temp.determinant();
}

/** Calcule la comatrice.
 *
 * La comatrice est la matrice des cofacteurs de l'instance courante. Chaque
 * élément \f$\left( i, j \right)\f$ est le cofacteur \f$\left( i, j \right)\f$
 * de l'instance courante.
 *
 * @return la matrice des cofacteurs de l'instance courante.
 */
Matrice3x3 Matrice3x3::comatrice() const
{
    Matrice3x3 temp;

    unsigned i, j;

    for ( i = 0; i < 3; i++ )
    {
        for ( j = 0; j < 3; j++ )
            temp( i, j ) = cofacteur( i, j );
    }

    return temp;
}

/** Opérateur unaire moins.
 *
 * Cet opérateur permet d'obtenir une nouvelle matrice \f$M^-\f$, où \f$ m^-_{i,
 * j} = - M_{i, j} \f$.
 *
 * Fait appel l'opérateur *= pour une nouvelle instance.
 *
 * @return une nouvelle matrice, égale à \f$ -1 \f$ fois l'instance courante.
 */
Matrice3x3 Matrice3x3::operator-() const
{
    Matrice3x3 temp( *this );

    return temp *= -1.;
}

/** Opérateur de transposition.
 *
 * Cet opérateur permet d'accéder à la matrice transposée \f$ M^T \f$, pour
 * laquelle \f$ M^T_{i, j} = M_{j, i}\f$.
 *
 * @return une nouvelle matrice, égale à la transposée de l'instance
 * courante.
 */
Matrice3x3 Matrice3x3::transpose() const
{
    Matrice3x3 temp;

    for (unsigned int i(0), j; i < 3; i++)
    {
        for (j = 0; j < 3; j++)
            temp( i, j ) = at( j, i );
    }

    return temp;
}

/**
 * Permet d'additionner une matrice à l'instance courante, qui est
 * modifiée.
 *
 * On délègue l'addition à la classe mère.
 *
 * @param[in] mat matrice que l'on veut additionner à l'instance
 * courante
 *
 * @return une référence sur l'instance courante.
 */
Matrice3x3& Matrice3x3::operator+=(const Matrice3x3& mat)
{
    matrix< double >::operator+=( mat );

    return *this;
}

/**
 * Permet de soustraire une matrice à l'instance courante, qui est
 * modifiée.
 *
 * On délègue la soustraction à la classe mère.
 *
 * @param[in] mat matrice que l'on veut soustraire à l'instance
 * courante
 *
 * @return une référence sur l'instance courante.
 */
Matrice3x3& Matrice3x3::operator-=(const Matrice3x3& mat)
{
    matrix< double >::operator-=( mat );

    return *this;
}

/** Permet d'effectuer la multiplication @e à @e droite par une autre instance.
 *
 * Cet opérateur permet de multiplier l'instance courante par une autre matrice
 * (multiplication @e à @e droite!). Le produit matriciel est défini comme :
 * \f[ \left( A \cdot B \right)_{i, j} = \sum_{k = 1}^3 A_{i, k} \, B_{k, j }
 * \f]
 *
 * On délègue en fait la multiplication à la classe mère en utilisant
 * la fonction prod.
 *
 * @param[in] mat matrice pour laquelle on veut faire la
 * multiplication @e à @e droite de l'instance courante.
 *
 * @return une référence sur l'instance courante.
 */
Matrice3x3& Matrice3x3::operator*=(const Matrice3x3& mat)
{
    //matrix< double >::operator*=( mat );
    *this = prod( *this, mat );

    return *this;
}

/** Opérateur d'addition.
 *
 * Permet d'obtenir le résultat de l'addition de deux matrices.
 *
 * On utilise la méthode operator+= pour une matrice temporaire.
 *
 * @param[in] mat matrice à additionner à l'instance courante.
 *
 * @return une nouvelle matrice, égàle à l'addition des deux matrices.
 */
Matrice3x3 Matrice3x3::operator+(const Matrice3x3& mat) const
{
    Matrice3x3 temp( *this );

    return temp += mat;
}

/** Opérateur de soustraction.
 *
 * Permet d'obtenir le résultat de la soustraction de deux matrices.
 *
 * On utilise la méthode operator-= pour une matrice temporaire.
 *
 * @param[in] mat matrice à soustraire à l'instance courante.
 *
 * @return une nouvelle matrice, égàle à la soustraction des deux
 * matrices.
 */
Matrice3x3 Matrice3x3::operator-(const Matrice3x3& mat) const
{
    Matrice3x3 temp( *this );

    return temp -= mat;
}

/** Opérateur de multiplication.
 *
 * Permet d'obtenir le résultat de la multiplication de deux matrices.
 *
 * On utilise la méthode operator*= pour une matrice temporaire.
 *
 * @param[in] mat matrice qui va multiplier l'instance courante par la
 * \e droite.
 *
 * @return une nouvelle matrice, égàle à la multiplication des deux
 * matrices.
 */
Matrice3x3 Matrice3x3::operator*(const Matrice3x3& mat) const
{
    Matrice3x3 temp( *this );

    return temp *= mat;
}

// Opérateurs sur des scalaires
/** Permet de multiplier la matrice par un nombre réel.
 *
 * On délègue l'opération à la classe mère.
 *
 * @param[in] scalaire nombre réel.
 *
 * @return une référence sur l'instance courante.
 */
Matrice3x3& Matrice3x3::operator*=( const double& scalaire )
{
    matrix< double >::operator*=( scalaire );

    return *this;
}

/** Opérateur de multiplication par un scalaire.
 *
 * Permet d'obtenir le résultat de la multiplication d'une matrice par
 * un nombre réel.
 *
 * On fait appel à la méthode operator*= pour une instance temporaraire.
 *
 * @param[in] scalaire nombre réel.
 *
 * @return une nouvelle matrice résultant de la multiplication de
 * l'instance courante par un nombre réel.
 */
Matrice3x3 Matrice3x3::operator*(const double& scalaire) const
{
    Matrice3x3 temp( *this );

    return temp *= scalaire;
}

Matrice3x3 operator*(const double& a, const Matrice3x3& mat)
{
    return mat * a;
}

// Opérateur sur des "vecteurs"
/** Opérateur de multiplication sur des points de l'espace.
 *
 * Permet de multiplier @e à @e doite l'instance courante par un point de
 * l'espace.
 *
 * Fait appel à la fonction prod.
 *
 * @param[in] pt point de l'espace (Coord) que l'on veut multiplier
 * par \e la \e droite par l'instance courante.
 *
 * @return un nouveau point de l'espace \f$ \vec{OP'} = M \, \vec{OP}
 * \f$.
 */
Coord Matrice3x3::operator*(const Coord& pt) const
{
    Coord res( prod( *this,  pt ) );

    return res;
}

/** Opérateur de multiplication sur des vecteurs de l'espace.
 *
 * Permet de multiplier à doite l'instance courante par un vecteur de
 * l'espace.
 *
 * Fait appel à la fonction prod.
 *
 * @param[in] vec vecteur de l'espace (Vecteur) que l'on veut
 * multiplier par \e la \e droite par l'instance courante.
 *
 * @return un nouveau vecteur de l'espace \f$ \vec{v'} = M \, \vec{v}
 * \f$.
 */
Vecteur Matrice3x3::operator*(const Vecteur& vec) const
{
    Vecteur pt( prod( *this, vec ) );

    return pt;
}

Coord& operator*=(Coord& pt, const Matrice3x3& mat)
{
    pt = prod( pt, mat );

    return pt;
}

Vecteur& operator*=(Vecteur& vec, const Matrice3x3& mat)
{
    vec = prod( vec, mat );

    return vec;
}

Coord operator*(const Coord& pt, const Matrice3x3& mat)
{
    Coord temp( pt );

    return temp *= mat;
}

Vecteur operator*(const Vecteur& vec, const Matrice3x3& mat)
{
    Coord temp( vec );

    return temp *= mat;
}

ostream& operator<<( ostream& os, const Matrice3x3& mat )
{
    for (unsigned int i(0), j; i < 3; i++)
    {
        for (j = 0; j < 3; j++)
            os << mat.at( i, j ) << ' ';
        os << endl;
    }

    return os;
}

/** Opérateur d'affectation.
 *
 * Opérateur d'affectation d'une instance de la classe mère vers
 * l'instance courante. Est défini pour simplifier le code interne.
 *
 * Le mécanisme de copie est laissé à la charge de la classe mère.
 *
 * @param[in] mat instance de Matrice3x3 à copier.
 *
 * @return une référence sur l'instance courante.
 */
Matrice3x3& Matrice3x3::operator=( const MatriceBase& mat )
{
    MatriceBase::operator=( mat );

    return *this;
}
